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The equations of relativistic hydrodynamics are transformed so that steps forward in time pre- 
serves local simultaneity. In these variables, the space-time coordinates of neighboring points on 
the mesh are simultaneous according to co-moving observers. Aside from the time step varying 
as a function of the location on the mesh, the local velocity gradient and the local density then 
evolve according to non-relativistic equations of motion. Analytic solutions are found for two one- 
dimensional cases with constant speed of sound. One solution has a Gaussian density profile when 
mapped into the new coordinates. That solution is analyzed for the effects of longitudinal acceler- 
• ation in relativistic heavy ion collisions at RHIC, especially in regards to two-particle correlation 

f"*) ' measurements of the longitudinal size. 
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£^ \ The large collective flow observed at the Relativistic Heavy Ion Collider (RHIC) have revitalized interest in rela- 
tivistic hydrodynamics for nuclear collisions [3, 0, 0, 0] • Numerical approaches for solving non-relativistic equations 
can be divided into numerous categories. Information regarding the densities and collective velocities might be stored 
for static mesh points (Eulerian) , or for mesh points that move with the fluid (Lagrangian) . Other approaches discard 
' the notion of a mesh and instead follow the motion of particles which represent fluid elements 0, [g, 0] , thus making 
the solution similar to that of molecular dynamics. Adding viscosity and heat conduction brings about an even wider 
range of approaches. For instance, once heat diffusion is introduced one must choose whether the velocity refers to 
that of the energy 0] (Landau frame) or that of the number density 0] (Eckart frame), or even of the entropy density 
[To| . In addition to specifying the character of the mesh and the equations, there exist innumerable algorithms. 
Relativistic modeling involves another class of choices, as there may be good reasons to consider meshes where 
. information is not stored for a space-time hyper-surface defined by a fixed Cartesian time, but rather for some 
""^ ' alternately defined hyper-surfaces. For instance, cosmological modeling might use the proper time, i.e., the time as 
measured by a local observer. For relativistic collision, where there is semi-Hubble-like flow along the z-axis, one 
1 1 ' might use the Hwo-Bjorken time [ill, T = Vt 2 — z 2 , which is the time as measured by an observer moving with 
constant velocity v = z/t. If the matter accelerates along the z axis, this time would differ from the time measured 
by an observer moving with the fluid. Once the fluid velocities significantly differ from that of the mesh, relativistic 
effects become important. 

In this study we investigate the prospects of propagating so that after each point i moves forward by 8ti, the 
hyper-surface will satisfy local simultaneity, 
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H ' u-dr = 0, (1) 

where dr^ is the differential distance between neighboring points on the hyper-surface, and u is the four-velocity, 
where in the frame of fluid u = (1, 0, 0, 0). Thus, in the frame of the fluid the time difference dt between neighboring 
mesh points is zero. After presenting the equations of motion that preserve local simultaneity in the next section, we 
present analytic solutions in these coordinates in the next two solutions. The first solution is for a Gaussian-shaped 
entropy distribution that is initially at rest, while the second is a restatement of the Khalatnikov-Landau solution for 
matter leaving a static shock [l3l Il4j . Both solutions are for the particular case of purely one-dimensional motion 
with a constant speed of sound. Sec. IIVI explores the manifestations of longitudinal acceleration and compares the 
analytic Gaussian case with a numerical solution for the case with initial collective flow given by the Hwa-Bjorken 
velocity profile, u z = z/t. In particular, the degree to which the velocity gradient differs from 1/r is discussed for its 
importance in affecting conclusions about lifetime from two-particle correlations. The prospects and limitations for 
these approaches are reviewed in Sec. IVI1 after which an appendix is included reviewing other analytic solutions to 
relativistic hydrodynamics. 



II. MAINTAINING SIMULTANEITY IN RELATIVISTIC HYDRODYNAMICS 



Here, we express the equations for relativistic hydrodynamics in terms of coordinates relevant for a set of observers 
traveling with the fluid. Each observer initializes his clock so that, in each observer's frame, evolution will begin 
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simultaneously for the neighboring observers. Such is the case for the commonly used Hwa-Bjorken assumption 
[lH where the collective velocity is given by u z — z/t, with t being the proper time, i.e., the time measured by 
each observer in their frame. In the Hwa-Bjorken picture, there is no acceleration and local simultaneity (u ■ x = 0) 
is preserved if each co-moving observer propagates forward equal amounts in proper time. However, as will be shown 
here, simultaneity will not be maintained in the presence of acceleration. For accelerating cases, each co-moving 
observer must propagate forward by different time steps. 

The equations of motion for a relativistic hydrodynamic fluid conserving entropy are: 

(u ■ d)s = — sd ■ u (2) 
(u-d)u a = -jj^-(d a -u a u-d)P, (3) 

where s is the entropy density as measured by a co-mover. Throughout this section, we ignore the general case where 
other currents are conserved, thus making P a function of s only. For a specific differential fluid clement, an observer 
moving in the frame of the fluid observes a change in the velocity and entropy density, 

5v = -^-Jr, (4) 
Ss = -s(V-v)5t (5) 

where St is the differential time step and r and v are the position and velocity measured by an observer moving with 
the same velocity as the matter. 

The equations of motion can be integrated forward numerically on a Lagrangian mesh, one that moves with the 
fluid, by solving for the time development of the separation between neighboring mesh points, Ar, and the relative 
velocity of the two points, Av, where Ar and Av are always measured by an observer traveling with the fluid between 
the two mesh points. Since the velocities are, by definition, small according to this observer, the motion appears non- 
relativistic. After solving for the non-relativistic evolution of Av and Ar, one could perform incremental relativistic 
velocity additions and Lorentz transformations to solve for the relativistic evolution. Thus, it would appear that 
relativistic complications could be avoided until after the non-relativistic evolution of the mesh is determined. 

The algorithm described above assumes that Ar represents the separation of two simultaneous space-time points, 
i.e., At = in the frame of the fluid. Unfortunately, even if the mesh is originally defined with neighboring points 
being simultaneous, the simultaneity will be destroyed after propagating each point forward a fixed amount of proper 
time if the points undergo acceleration. Here, we show that simultaneity can be maintained if different points on the 
hyper-surface are propagated forward by different steps in proper time. 

Simultaneity between two neighboring points x\ and X2 with velocities u\ and u 2 can be stated by the relation, 

(ui +ita) • (x2 - %x) = 0, (6) 

where Ax = X2 — x\ is small and (1/2) (u\ + u 2 ) is the four velocity of a co-moving observer between the two points 
of the fluid. Maintaining simultaneity requires: 

(ui + u 2 + clSti + aSr 2 ) ■ {xi — x 2 + u\5ti — u 2 5t 2 ) = 0, (7) 

where a M is the acceleration. Since accelerations have no time components in the frame of the matter, a • u = Q, and 
since (u\ — u 2 ) ■ {u\ + u 2 ) is also zero, one can rewrite Eq. ([7]) to find: 

[(ui + u 2 ) + a(Sri + St 2 )} ■ [Ax + {m - u 2 )(Sn + St 2 )/2 - (m + u 2 ){8t 2 - 5n)/2] = 0, (8) 

(S Tl + 5r 2 )a ■ (Ax) - 2(St 2 - 5 Tl ) = 0, 
8t 2 = <5ti(1 — a ■ Ax), 

where it has also been assumed that Ax and St are small, i.e., (ui + u 2 ) ■ (u\ + u 2 ) « 4. The last expression in Eq.s 
© can be integrated to find an expression for ratio of time steps between an observer at x and a reference observer 
at x R , 

(9) 
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where the last step used the definition of the speed of sound for the case where P is expressed as a function of the 
entropy density, 

Z - 4*7 do) 

Thus, although the equations of motion, Eq.s (f?)), appear simple to solve by working in a co-moving frame where 
distance is measured along a path satisfying u-dx = 0, the equations can not be integrated forward by fixed St without 
violating the condition of simultaneity. Instead, one should integrate forward in time by an amount, dr = a(x)dr^, 
where tr refers to the time measured by a reference observer at xr. 

The equations of motion can be rewritten in terms of dru rather than dr, 

5s = -as(V • v)5t r , (11) 
P + e 

These expressions differ from the non-relativistic equations of motion due to the presence of a on the right-hand sides 
of both equations. These equations are especially well suited for Eulerian approaches, where the mesh moves with the 
fluid, and one need only solve the equations of motion for the evolution of differences of position and velocity between 
neighboring points, as measured by a co-moving observer. 

Solution in more than one dimension of Eq. (|11[) requires that the flow is ir-rotational. For rotational flow, the 
definition of simultaneity, u ■ dx = 0, can not uniquely map to a hyper-surface extending over all space. For instance, 
one can consider a ring of particles in uniform circular motion of a radius R at a velocity v. If two neighboring 
particles would synchronize their clocks as to be simultaneous in the frame of an observer moving with the particles, 
the laboratory observer would view this synchronization event as happening at two different times separated by 
St = vAx. If the synchronization were to be performed clockwise around the ring, the times would be off by a factor 
2ttRv after returning to the original point. The problem with rotational flow is also encountered when considering for 
a = exp(— J dx ■ a) which is only path independent if V x a = 0. If V x v = everywhere for the initial condition, and 
if V x a — for all times, the flow will remain ir-rotational. The requirement that V x a — is maintained when the 
pressure depends solely on the energy density, but can be violated if the pressure also depends on a particle density 
7i. This can be shown by considering 

Vx* = Vx-^L (12) 
P(e, n) + e 

Ve x VP 

~ {P + e? ( ' 
dPVexVn 

~ ~~8t~1 {P + e) 2 ' (1 ' 

For the physics of RHIC, and especially for the physics of the LHC, the particle densities are small and rotational 
aspects of initial flow should be small. The flow would indeed be rotational at non-zero impact parameters, even 
at midrapidity However, these effects should be less pronounced at higher beam energies if the production 

mechanism becomes increasingly boost-invariant. 

III. ONE-DIMENSIONAL GAUSSIAN SOLUTION 

Here, we present an exact analytic solution to the hydrodynamic equations of motion in Eq. (jlip for the specific case 
of one-dimensional motion where the speed of sound is constant. The solution is found by first mapping the coordinate 
£, which measures the distance along a space-time path of simultaneity, to a new coordinate I. The transformation 
to the new coordinate I will be chosen to cancel the factor of a in the equations of motion. 



d£ = ad£, £= I d£'a{£'). (15) 
In these coordinates, the equations of motion become 



8s = —s — s Str, (16) 
d£ 

Sv = — Str. 

P + e d£ 
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Expanding 5s and 5v, 

r)e <9q r\9 fin 

(17) 
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Here, we have identified the sum of small changes of the velocity as measured by co-movers, Ai>, as the rapidity y. 
These would look exactly like the non-relativistic expressions if the velocity could be identified with d£/dr R . This will 
not always be satisfied. In general, one can write £ as a sum over contributions between neighboring cells of length 
At, 

£ - ^ !>«»«<) < 18 » 

i 

EdAi y-^ da % 



^— ' dr ^— ' dr R 

f ^ da t 

The last term will vanish only if the co-moving time derivate of a vanishes. This criteria will be satisfied for our 
specific case where the speed of sound is constant and the entropy profile is set by a scaling function. From the 
definition of a in Eq. (J9j) one can see that for a constant speed of sound, 

a = exp {c 2 s Hs/s R )} = f-^j ' . (19) 

For a scaling solution, s(£, tr) = F(1/R(tr)), and the shape of the entropy profile is unchanged through time, aside 
from changing the scale R. The ratio s/sr will be then stay constant for each co- mover. 

Here, we present a solution where the entropy profile is assumed to follow a Gaussian profile, 

^"■'-^""{-ssra 1 }- (20) 

Furthermore, the velocity profile is assumed to be linear in £, 

?L = A{t = c s T R /Ro)^-. (21) 
c s Ro 

Since the entropy profile scales with time, the equations of motion are those of non-relativistic hydrodynamics, 

ds d(sy) 



dr R dl 



(22) 



dy dy 2 d In s 

d^ V ~d£ Cs dl 

After applying the hydrodynamic equations of motion to the assumed forms for s and y, the dependence on £ factors 
away and leaves simple differential equations of motion for A(t = c s t r /Rq) and for r(t) = R(t)/Ro, 



dA .o 1 
at r z 

dr/dt 



A = 

Substituting for A provides a differential equation for r, 

d 2 r 1 



dP r' (24) 
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which can be solved by making an analogy with classical mechanics where one wishes to find the trajectory for a force 
behaving as 1/r. The solution is 



t = I _^= = 2 1 /2 / ■ dy e y\ (25) 




1 y/2\li{x) 

where D(x) is Dawson's integral function [l6| defined by the integral, 

D(x) = e~ x ' 2 f du e u \ (26) 
Jo 

D(x) rises proportional to x for low x and falls as l/2x for large x. Here, t = refers to the time at which the velocity 
gradient vanishes, and r represents the ratio of the R to the radius at that time. The expression is analytic, but in 
practice, evaluating the expression might require programming the function, which diminishes the glamour of finding 
an analytic expression. This expression must be inverted to find r as a function of t. The velocity gradient can then 
be expressed in terms of r, 

1 



A (i) = -r-VH^) (27) 
r(i) 

A disquieting aspect of the solution is that the entropy density is s = dS/d£, while the solution conserves the 
integral J sd£. However, one can understand that both dS/d£ and dS/d£ represent conserved charges by considering 
a co- moving cell. The two quantities differ by a factor a — d£/d£, and since a remains constant and since the entropy 
is conserved within each cell individually, both J sd£ and J sd£ survive as constants of the motion. 

If one solves the equations of motion in terms of tr and £, the state of the system in terms of the coordinates 
determined in a fixed reference frame, x and t, requires inversion of the mapping procedure. To find s and u for a 
given time and position in the x direction, one could follow the following steps: 

(i) As a function of £ and tr, use the equations of motion to solve for the entropy density, s, the rapidity, v, and the 
factor a. 

(ii) Begin with tr — t and £ = 

(iii) For small Aa;, choose Ar and A£ to satisfy the conditions, 

Ax = cosh(y)A^ + sinh(y)Ar, (28) 
= cosh(y)Ar + sinh(y)A£. 

(iv) Find A£ = aA£ and Atr = oAr. Increment tr and £. 

(v) Repeat steps (iii) - (iv) until the point (t, x) is reached. 

Solutions for different values of i?o and c s can be transformed into one another by scaling tr and R through the 
dimensionless variables t and r described above. However, this simple equivalence is lost when one maps to the 
coordinates in fixed reference frame, x and t. The mapping procedure will then depend on the ratio of the speed of 
sound to the speed of light. 

Figure [1] shows the entropy density and collective rapidity as a function of the position and time in a Cartesian 
reference frame for the case where the speed of sound is c? s = 0.2. In terms of £, the profile is of a Gaussian form, 
but after scaling to d£ = di/a, then performing the Lorentz transformations to Cartesian space-time coordinates, the 
profile is decidedly non-Gaussian. This is evident for the entropy density, s, plotted in the upper panel of Fig. [T]for 
t = 0. On a logarithmic scale, a Gaussian is concave downward for all x; whereas the scaled Gaussian picks up a 
power-law tail, which after some non-trivial expansions can be shown to have the form, 

/ \-l/c 2 

s(x -> oo, t = 0) ps ( — == ) , u = c s x/{R V2). (29) 



Since the acceleration for a fixed speed of sound is proportional to the logarithmic derivative of the density, the 
reversal of the concavity in s(x) results in an acceleration that is maximized for intermediate values as illustrated by 
the plot of rapidity vs. position for subsequent times displayed in the lower panel of Fig. [TJ This is distinctly unlike 
the behavior one would observe for a non-relativistic system [l7^ . where the acceleration rises linearly at all x. 
The solutions illustrated in Fig. [T]were also checked for conservation of energy and entropy. 

E t ot(t) = [ dxT ao (x,t), S tot = / dx u°{x,t)s(x,t), (30) 
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FIG. 1: The analytic solution is shown after being transformed into Cartesian space-time coordinates for three times, t — 
(circles), t = 10-Ro/cs (squares) and t — 20Ro/c s (triangles). The initial entropy profile is Gaussian in the variable £ but in 
Cartesian coordinates behaves like a power-law in the tail. The velocity profiles turn over at large x is due to the falling slope 
of the entropy profile combined with the effects of transforming the solution into Cartesian coordinates. 



where the stress-energy tensor was calculated from s(x, t) and the collective four velocity w, 

T al3 (x,t) = {P + ()u a uP - Pg af3 , (31) 
P = s 1+c % e = P/c 2 s . 

Although the expressions used to calculate u and s are analytic, calculation of the Dawson's integrals, and the mapping 
to a fixed reference frame involved numerical integrations which introduced small errors. These errors were less than 
0.1%, and fell with increasing granularity in the calculations. Furthermore, the net momentum of all matter moving 
to the right was tested for consistency with the integrated impulse, 



J dx T 0x {x,t) = j dt'P{x = 0,t'). 



(32) 



This condition was also satisfied to high precision. 

Despite the fact that the space-time evolution of the entropy density is rather complicated, the entropy per rapidity 
has a simple form. To see this, we begin with the fact that s = dS/d£, then apply the chain rule to find dS/dy, 

dS__ Mdl _ 1 [ g \ dl d£ 

dy~ S Mdy~ R{T R ) e * V \ 2R{ TR y\ didy' ( } 

After using the simple linear form, y = AI/Rq, in Eq. (|2~lj) , and integrating Eq. ([9]) to find, 

a = d£/d£ = e- c * p /V R2 \ (34) 
After substituting for £ using Eq. (f2Tj) the rapidity distribution becomes 



dy V^F) I 2A2c l r2 



1 -pI^ITs m 
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Here, we have also substituted for rA = [ln(r 2 )] 1//2 using Eq. (|2T[) . and r is again related to R by scaling, r = R/Rq. 
The Gaussian width of the rapidty distribution is thus, 



collective velocity, R\ ollg dv/dz, surpasses the thermal velocity 



ln(r 2 ) f 

^ = Cs VT^f (37) 

For large times, the width of the rapidity distribution in Eq. ([35]) , grows logarithmically with time as can be seen by 
considering the large r behavior of the Dawson's integral function in Eq. (j25|) . y/ln(r) ~ ^/ln(t). We emphasize that 
this simple Gaussian form relies on the assumption that the breakup is sudden along a hyper-surface of simultaneity. 

Experimentally, the breakup time is inferred from two-particle correlations through the measurement of i?i on g which 
represents the size along the beam axis of the phase space distribution of particles with zero rapidity. If the collective 
motion covers large amounts of rapidity, the size of the phase space region is defined by the point at which the 

iiy 0.[hJ. 

p ^thcrm f na\ 

Rlons ~ d^Jd~z- (38) 

For the Hwa-Bjorken solution, where there is no longitudinal acceleration, the collective velocity is always z/t. The 
velocity gradient is then 1/r, where r is the proper time elapsed since the beginning of the collision. In the no- 
acceleration Hwa-Bjorken ansatz, the time r is then inferred by 

r = ^L. (39) 

^thcrm 

Here, i?i n g is experimentally determined and the thermal velocity is inferred from spectra or blast-wave models. The 
simple statement above is modified by collective transverse flow, though the basic linear nature of the dependencies 
is preserved. 

The presence of longitudinal acceleration would alter the inference of r by voiding the Hwa-Bjorken equivalence 
dv/dz = 1/t. In the analytic Gaussian solution shown here, the velocity gradient is given by Eq. (f2"T)) , 

dv 



dz ^V2HRm (40) 



1 



= -[2xD(x)\, x = a yy /(l - c 2 )/2c 2 . 

T 

Here, the second line used the relation for the rapidity width, Eq. (j37|) , and the final step used Eq. ([25]) . 

The quantity in the square brackets in Eq. (|40|) now represents the correction to the Bjorken formula for the 
lifetime, 

t = \2xD{x)] . (41) 

^thcrm 

For large rapidty widths, the Dawson's integral function can be expanded for large x, 

1 1 3 (271-1)!! 

The first order term reproduces the Hwa-Bjorken result. At RHIC, the variance of the pion rapidity distribution 
is a bit over 2.0 (20j . However, here y refers to the source rapidity which should combine in quadrature with the 
thermal spread, ~ 1 to result in the pion rapidity spread. Thus, a v should be ~ 1.7 here, which if we combine this 
value of a with a speed of sound, c 2 ~ 0.20, x should in the neighborhood of 2.4, which would suggest that the 
Hwa-Bjorken estimate underestimates the lifetime by approximately 10%, which is similar in strength to what was 
seen by comparing three-dimensional to boost-invariant numerical hydrodynamics by Renk, [2l| . 



IV. COMPARISON WITH HWA-BJORKEN INITIAL CONDITION 



The analytic solution of the previous section was predicated on the physical picture of complete stopping. Ex- 
periment strongly suggests, mainly through baryon distributions, that there is significant stopping, which invalidates 



8 



the analytic expressions derived in the previous section. However, it is often the case that for large times solutions 
approach analytic forms even though the initial conditions differed significantly from the analytic form. To investigate 
this possibility, we compare the solutions for the same speed of sound, c 2 s = 0.2, but for two very different initial 
conditions: 

(a) An analytic solution which assumes complete stopping at t = 0. An initial Gaussian size of 0.06416 fm was 
chosen so that when the system expanded to its breakup time of 10 fm/c, the rapidity width would match that 
estimated for the final state as discussed in the previous section, a y (r = lOfm/c) = 1.7 [26j. Here, the variance 
is defined using the entropy density as a weight, 



fd£y 2 s(£,T) 



(b) Numerical solutions of hydrodynamics assuming a Hwa-Bjorken ansatz for the initial conditions, with the ra- 
pidity at time To = 0.5 fm/c given by, 

y = (44) 

TO 

where £ is defined in the same way as in Sec. [TT] as the net distance separating co-movers. The initial entropy 
density is chosen to be Gaussian, 



dS_ 
~d£ 



oc e 



(45) 



These initial conditions differ from that used in the analytic solution both in that the matter is not moving, and 
in that the Gaussian form was applied to £, rather than the scaled variable £. The initial time To was chosen 
rather arbitrarily to be 0.5 fm/c, and the initial size i?o was then adjusted so that the variance a y would again 
be 1.7 at a breakup time r = 10 fm/c. The size required to achieve this fit is Rq = 0.55 fm. 

Thus, the two solutions are matched by having the two solutions break up at the same reference time along a 
hyper-surface of simultaneity, and by having the same rapidity spread in the final state a y . Figure [2] demonstrates 
the remarkable agreement of the two evolutions for all but the earliest times by showing a y and the inverse velocity 
gradient at the origin as a function of r. At r = 0, the analytic solution has zero velocity gradient, but by the time the 
other solution begins at t = 0.5 fm/c has already reached a similar amount of collective velocity. Figure [3] shows that 
the entire structure of the density and rapidity profiles are very similar at breakup. This shows that the final state 
of the system tends to move towards the analytic solution at large times, with a given final state being reached by a 
large number of possible initial conditions. We would expect this trend to be violated for when comparing solutions 
with very different equations of state. 

Although these models have an overly simplistic equation of state and ignore the effects of transverse flow, one can 
gain some quantitative understanding of the effects of longitudinal expansion at RH1C from the evolutions shown 
in Fig. [21 First, the variance of the collective rapidity a y grew by more than 60% between r = 0.5 fm/c and 10 
fm/c. Even after incorporating transverse flow or using a softer equation of state, the growth would be several tens of 
percent. This clearly plays a critical role in understanding the stopping at RHIC, as ignoring the effects of longitudinal 
acceleration would lead to an underestimate of the initial stopping in central collisions. 

The evolution of the inverse velocity gradient at the origin in the lower panel of Fig. [5] provides quantitative insight 
into how longitudinal acceleration affects the determination of the lifetime from correlations. In correlation analyses, 
one obtains the inverse velocity gradient by dividing the longitudinal correlation radius, i?i on g, by the thermal velocity 
^therm, which is determined from analyses of spectra. Thus, the analyses first determine the inverse velocity gradient, 
which are identified with the breakup time by assuming an acceleration-less Hwa-Bjorken form for y vs. £. Fig. [2] 
shows that the neglect of acceleration leads to an underestimate of the break-up time of ~ 10%. As discussed in 
the previous section, Eq. (|4ip shows that this underestimate depends only on the ratio of the speed of sound to a y 
at breakup. Given that solutions tend towards the analytic solution, this conclusion is somewhat robust. However, 
lowering the speed of sound or incorporating transverse flow should ameliorate this effect. 

Some analyses of correlation and spectra have been based on assumptions of a purely Boost invariant picture for 
the longitudinal flow. This leads to an underestimate of the lifetime for two separate reasons. First, the system is not 
infinite in extent. Convoluting a boost-invariant solution with a Gaussian profile should lead to larger estimates of the 
lifetime, with the amount of the correction being of the order of 10%, depending on the p f range studied. The second 
effect is the one discussed here, of longitudinal acceleration. Boost-invarian blast waves tend to suggest breakup times 
near 9 fm/c [22l |. These times give one pause since the same analyses suggest that the radii have grown from 6 fm 
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initially to 12 or 13 fm, with final collective velocities of 0.7c at the edge. To achieve this final state geometry, the 
matter would have had to achieve maximum transverse speed immediately, which is unphysical. This inconsistency, 
which represents part of the "HBT puzzle" , would be softened by increasing the estimates of the breakup lifetime to 
11 or 12 fm/c. 




r (fm/c) 

FIG. 2: The rapidity width and inverse velocity gradient are shown as a function of r for the analytic solution (solid lines) and a 
numerical solution of a system whose initial collective flow at to = 0.5 fm/c was that of a Hwa-Bjorken model. For the analytic 
solution, the matter was initially stopped with a size chosen so that the rapidity width would grow to 1.7 by 10 fm/c. For 
the numerical solution, the initial width of the rapidity distribution was chosen so that it would also reach a y = 1.7 at r = 10 
fm/c. Even though the initial velocity gradient is zero in the analytic solution, the models behave similarly for all r > r . The 
dashed lines show what would be expected in a Hwa-Bjorken expansion that neglects acceleration. The upper panel shows 
that acceleration increases the rapidity width by over 50%. In the lower panel, the difference between the calculations with 
acceleration and the acceleration-less solution demonstrates the error associated with assigning the inverse velocity gradient, 
which is inferred from correlation measurements, with the lifetime. Accounting for acceleration should increase the estimate of 
the breakup time by ~ 10%. 



V. ONE-DIMENSIONAL EXPONENTIAL SOLUTION 

Like the first part of the Khalatnikov-Landau solution, we will consider matter ejected from a slab of uniform energy 
density, thus re-deriving their result but using the mapping procedure applied to the Gaussian profile in the previous 
section. Unlike the Gaussian solution, all the matter undergoes acceleration in this solution, including the reference 
point. We assume a simple form for the rapidity profile and entropy density: 



y(£, t r ) = y R (r R ) + A(r R )£, s(I, r R ) = SR (r R )e-^ R( ^ . 



(46) 
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FIG. 3: The rapidity and entropy profile are shown as a function of the distance along a hyper-surface of simultaneity. Distances 
are scaled by the time r for easy comparison with a Hwa-Bjorken expasion. The analytic solution (solid lines), which assumed 
complete stopping at r = 0, is compared to a numerical solution for Hwa-Bjorken initial conditions with to = 0.5 fm/c. By 
t — 10 fm/c (squares), the two solutions are in almost perfect agreement. Even at r = 0.5 fm/c (circles), the two solutions 
match well for I < 3r. 



Here, tr is the time as measured by a reference observer who moves along with the fluid whose local entropy density 
is sr and moves with rapidity i/r. The distance £ is the scaled distance measured relative to the reference observer, 



s 



£ = J d£'a(f), a(£)=expj-y d£aj = ^—j , (47) 

where £ is measured along a line of simultaneity, u ■ dx = 0. Here, a is again the acceleration of the fluid element. 
This parameterization will solve the hydrodynamic equations of motion if 

1 B 

A{tr) = — , R(t r ) = c s t r , sr(tr) = — (48) 

TR TR 

where c s is the speed of sound, and B is an arbitrary constant. As was the case with the Gaussian solution a stays 
fixed for a co-mover in this solution which allows one to state 

-^ = V-yn. (49) 

Since these solutions are in terms of positions and times relative to an accelerating observer, expressing the solution 
in Cartesian space-time coordinates requires first solving for the trajectory of the reference observer. This is not 
difficult as the observer has velocity c s at t = and accelerates with a = 1/tr after that point. Since these solutions 
have the same initial conditions and satisfy the same underlying hydrodynamic equations of motion, they should 
reproduce the first-stage Khalatnikov-Landau solution. This was checked numerically. 

The rapidity distribution can be extracted from the solution by considering dS/dy which should approximately be 
proportional to the number of particles per unit rapidity. This is only approximate because thermal motion smears 
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out the correspondence between a particles asymptotic rapidity and the rapidity of the source from which it was 
emitted. From the solution above, 

dS = g d£d£ = Be _ {y _ yR){1/Cs _ Cs) 
dy dydl 

where dijdl = a = e - ci / TR . The precise form of dS/dy depends on the form of the breakup surface in space-time, 
which is defined by the simultaneity condition here. A remarkable aspect of this result is that the shape is determined 
solely by the speed of sound, and is independent of the breakup time. Of course, after the shocks meet, matter is no 
longer being fed into the solution at y = 0, and the width of the rapidity distribution would no longer be fixed. 



VI. CONCLUSIONS 



Here, we have investigated the prospects of solving hydrodynamics in a coordinate system that maintains local 
simultaneity. This involves storing information about the number and entropy densities and about the velocity gradient 
along a hyper-surface defined by the condition u-dx = 0. At each point in time the subsequent evolution of the velocity 
differences to the next time step become equivalent to those of non-relativistic hydrodynamics. However, maintaining 
simultaneity requires propagating by a time step St that differs at each mesh point by a factor a = exp(— J a ■ dx) 
relative to a reference observer. 

This approach has both advantanges and limitations. One advantage is that the differential equations are linear in 
the velocity, unlike the usual equations which involve 7 factors. Since mesh points satisfy u ■ dx = 0, it is particularly 
convenient to extract information about the matter in its rest frame from the mesh. The coordinate system is especially 
attractive for one-dimensional systems, where for certain cases, one can also scale the position measurement and find 
some exact analystic solutions by mapping to non-relativistic solutions. 

Working with Simultaneity-conserving coordinates does involve added difficulty when compared to non-relativistic 
hydrodynamics. First, one must calculate the coefficient a at each mesh point. For ideal hydrodynamics this is not 
difficult as the coefficient depends only on the local densities, and the densities of the reference observer. For viscous 
hydrodynamics, calculation of a might involve performing a path integral of the acceleration. Viewing the evolution 
in standard coordinate systems, e.g., a Cartesian space-time mesh, can also be clumsy as the transformation to a 
different hyper-surface is a bit cumbersome. 

The most serious limitation of the method concerns rotational flow. In the absence of rotational flow the condition 
u • dx = defines an infinite hyper-surface in space time, and each step forward in time defines a new surface. After 
evolving forward in time, each point in space time will belong to one surface, and one surface only. This is not the case 
for rotational flow, as the condition of simultaneity becomes path dependent. This restriction certainly disqualifies 
this approach for many applications for relativistic heavy ion collisions. 

For certain one-dimensional cases, the method leads to analytic solutions. These are scaling solution, i.e., ones 
where the density profile is determined by a single parameter. In these cases the hydrodynamic equations of motion 
exactly reproduce the non-relativistic ones after the distances along the hyper-surface are scaled along with the time. 
Along with the Khalatnikov-Landau solution for matter ejected from a shock, which was already well known, we were 
able to find a Gaussian solution. In the scaled coordinates, it involves the expansion of a Gaussian-shaped profile, 
exactly like the well-known non-relatistic case. One of our biggest disappointments is that the analytic solutions 
were confined to purely one-dimensional solutions. We could not even find analytic solutions for problems with 
spherical symmetry. Analytic solutions in higher dimension would have been significantly more useful for testing 
highcr-dimcnsional relativistic codes. 

One virtue of the anlaytic solution was that it provided a transparent means for understanding the effects of 
longitudinal acceleration at RHIC. The analytic solution, which assumes complete stopping, and a numerical solution 
of a system with an initial Hwa-Bjorken expansion were compared with the constraint that they acheived the same a y 
at breakup. The analytic solution had such significant acceleration, that by a time of 0.5 fm/c, when the numerical 
solution was initialized, it had already achieved a similar amount of collective velocity, and by breakup the two 
solutions were indistinguishable. For a speed of sound c 2 s = 0.2, the rapidity spread was shown to increase by over 
50% between r = 0.5 fm/c and a break-up time of 10 fm/c. The quickness with which the initial conditions became 
irrelevant and the numerical solution merged with the analytic one was remarkable, and represents an important 
lesson. 

The effect of longitudinal acceleration is especially important for determining the lifetime from correlation analyses. 
The effect of the the expansion was to increase the life time by a simple factor, 2x/D(x), where x depended only on 
the ratio of a y at breakup to the speed of sound. For RHIC this suggests a 10% effect. Since the speed of sound and 
final rapidity spread are well constrained by other measurements, the correction factor factor is quite robust. The 
main qualifier comes from the fact that the effects of transverse expansion were not coupled into the evolution. 
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After having investigated prospects in detail, we only enthusiastically endorse simultaneity-conserving approaches 
for one-dimensional problems, or for problems for which the symmetry reduces the problem to one dimension. We do 
not foresee major technical issues for extending the approach to more complicated equations of state or for incorpo- 
rating viscous effects. For higher dimensions, the difficulties associated with the path integrals and the limitations of 
ir-rotational flow might bring along more trouble than is warranted by the advantages of the approach. 



Appendix: Other Analytic Solutions 

Here, we review three solutions to one-dimensional relativistic hydrodynamics. The first is the Khalatnikov-Landau 
solution [13l Il4j which assumes initial conditions of a thin slab with a uniform energy density inside and zero energy 
density outside. Since no information can travel into the slab, due to the fact that d x P — 0, a shock develops. In the 
frame of the discontinuity the matter moves outward from the front at the speed of sound c s . Choosing a coordinate 
system where the front is at x = 0, the velocity v and entropy density s outside the slab are described by: 

v{x,t) = £±£^, s(x,t) = aoe-y&M*', (51) 
t + c s x 

where y = tanh^ 1 ^) is the rapidity and So is the entropy density at x — which is determined by solving the 
Rankine-Hugoniot equations and is determined by the equations of state and the energy density inside the slab. This 
solution requires the additional assumption that the speed of sound is fixed. Once the two shock fronts meet, the 
solution becomes a real tour de force, and has a much more complicated form than the simple expressions above. The 
latter part of the solution is expressed in terms of derivatives of an integral of Bessel functions, which can be easily 
evaluated numerically. 

The second solution is the boost-invariant solution of Hwa and Bjorken [ll], [l^] is especially simple and can be 
easily incorporate an arbitrary equation of state. Contrary to the Landau parameterization, this solution assumes the 
existence of an initial flow velocity, and in fact, assumes that the flow covers an infinite rapidity window. The solution 
assumes that all frames which pass through the z = t = space-time point will see the same evolution. Referred to 
as "boost invariance" , this implies that v = z/t, and implies zero acceleration. Furthermore, the energy density e, the 
entropy density s, the temperature T and the pressure P depend only on the time measured by co- moving observers, 
t = 7< — jvz = \/t 2 — z 2 , which is known as the proper time. In the Hwa-Bjorken solution, the distance between 
neighboring points at the same proper time, as measured by a co-moving observer, is Sz — rSrj, where 

z = r sinh(f?), t = rcosh(?7). (52) 

Given that v = z/t, r) equals the rapidty y in the Hwa-Bjorken solution. The collective velocity gradient is the same 
for all co- moving observers, dv/dx — 1/r. 

Current conservation for the boost-invariant solution, d^j^ = 0, governs the time evolution of any conserved charge 
density. Since V • v = 1/r, current conservation states that all charge densities fall inversely with time, 

P = Po — ■ (53) 

T 

Since entropy is also conserved, s also falls inversely with r. Since P and e are functions of s and the charge densities, 
their time dependence is determined by the time equation of state. For instance, if the speed of sound is constant, 
P = P (s/s ) 1+c ° , and P = Po(t /t) 1+c = . The boost-invariant nature of the solution can easily incorporate arbitrary 
equations of state, and can also incorporate viscosity. When viscosity is added, entropy is no longer conserved but 
one can solve for s as a function of t using the Navier-Stokes equation [23l |24|] . 

The Hwa-Bjorken solution would seem justified in the case of extremely high beam energy, where the source covers 
many units of rapidity. In that limit, the pressure depends only slowly on 77, and since accelerations are driven by 
pressure gradients, the velocities of particular regions of the matter are fixed, thus justifying the condition v — z/t. 
Boost invariance also seems natural from the perspective of having matter created by the classical fields originating 
from the ultrarelativistic charges belonging to the incoming hadrons. For an observer at mid-rapidity, the currents 
appear as qc or —qc, where c is the velocity of light. If the observer boosts to a new frame, the currents will still 
appear to be qc and —qc as long as the the observer does not approach the speed of the colliding hadrons. If the 
interactions are determined by currents, one would expect an approximate boost invariance at mid-rapidity. 

The shortcoming of the Hwa-Bjorken solution comes from ignoring longitudinal acceleration. As we will see in 
the next section, this shortcoming could result in 10-20% changes in the interpretation of some observables. One 
interesting aspect of boost-invariant solutions is that the conservation of total energy is not enforced. Each fluid 
element expands and does work PdV . However, that work never goes into increasing the collective kinetic energy of 
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the fluid elements, since they all coast without acceleration. Instead, the central cell does work on the adjacent cell, 
which then does even more work on the next cell. For an observer at rest, the work done onto subsequent cells at time 
t increases exponentially, in a kind of covariant procrastination, until eventually, one reaches the light cone. Again, 
as long as many units of rapidity are covered, and as long as observations are confined to mid-rapidity, the overall 
conservation of energy can be ignored. 

Recently, Csorgo et al. |25[ published a new set of exact one-dimensional analytic solutions for hydrodynamics. 
Unlike the Hwa-Bjorken solutions these solutions incorporate longitudinal acceleration. For the Csorgo solution the 
coordinate i] = tan -1 (z/£) no longer equals the rapidity. Instead, 

V = Ary, (54) 

with A being independent of 77 and r. A second simplyfing assumption made in [25j is that P depends only on r. 
Using Eq. (|54p one can express the collective velocities and derivatives in the frame of an observer at 770 moving with 
V = Vo, 

d^ = (d T ,(l/r)d v ), (55) 
u M = [cosh(Af7 — 770), sinh(Ar? — 770)]. 

Here, the square brackets are used to represent the two relevant components of the four vector. The divergence of the 
velocity is then, 

d-u=-. (56) 

T 

Entropy conservation, dE = —PdV, can equivalently be stated in continuous variables as d T e = — (P + e)(d-u), which 
gives 

^=-^P + e). (57) 

The equations of motion for hydrodynamics, 

(M'9K = -p^(^-ti"(M.9))P, (58) 

after insertion of the expression for u M , and the assumption that P depends only on r become, 
A dPIdr dPIdr 

- [sinh 2 (777) , sinner/) cosh^)] = — - — [1,0]- — - — [cosh 2 (777) , sinh^Ti) cosher;)] , 7 = (A - 1). (59) 
t P + e P + e 

Both components of the equation will be satisfied either if A = 1, or 

ap \ 

(P + e). (60) 
Dividing this by the similar expression for the energy density, Eq. (|57[) . gives the speed of sound 

= 1- (61) 



o dP 



fixed S 



Thus, the Csorgo solution is rather particular, as it only applies when the speed of sound equals that of light. If the 
dimensionality is d, the speed of sound must then be c 2 = 1/d [27j • 
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